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Abstract. This is a computational study of gravity-driven fingering in- 
stabilities in unsaturated porous media. The governing equations and cor- 
responding numerical scheme are based on the work of Nieber et al. [Ch. 23 
in Soil Water Repellency, eds. C. J. Ritsema and L. W. Dekker, Elsevier, 2003] 
in which non-monotonic saturation profiles are obtained by supplementing 
the Richards equation with a non-equilibrium capillary pressure-saturation 
relationship, as well as including hysteretic effects. The first part of the study 
takes an extensive look at the sensitivity of the finger solutions to certain 
key parameters in the model such as capillary shape parameter, initial sat- 
uration, and capillary relaxation coefficient. The second part is a compar- 
ison to published experimental results that demonstrates the ability of the 
model to capture realistic fingering behaviour. 
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1. Introduction 

The transport of water and dissolved contaminants witliin the vadose zone is extremely 
important in a wide range of natural and industrial applications including protection 
of groundwater aquifers, irrigation, flood control, and bioremediation, to name just a 
few. Many of these applications exhibit preferential flow in which gravitational, viscous 
or other forces initiate instabilities that propagate as coherent finger-like structures. In 
fingered flow, water is able to bypass a signiflcant portion of the porous matrix and thus 
penetrate more rapidly than would otherwise be possible for a uniform wetting front; as a 
result, flngering can have a major impact on the transport of contaminants carried by an 
inflltrating fluid. A clear understanding of fingering phenomena can therefore be essential 
in the study of certain applications such as groundwater contamination. 

We focus in this work on preferential flow that is driven by gravitational forces arising 
from the difference in density between invading water and displaced air. The structure of a 
typical flnger consists of a nearly saturated "tip" at the leading edge, behind which follows 
a "tail" region having a uniform and relatively low saturation (see Fig. 1). As the flnger 
penetrates into the soil, the region immediately behind the tip drains somewhat causing 
pressure to decrease and preventing the finger core from widening, thereby allowing the 
unstable finger to persist in time. Experimental studies have provided additional insight 
into the detailed character of fingers and the physical mechanisms driving their formation, 
beginning with the work of Hill and Parlange [1972] and continuing to the present day 
with the work of authors such as Diment and Watson [1985], Glass et al. [1990], Selker 
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et al. [1992a], Lu et al. [1994], Banters et al. [2000], Yao and Hendnckx [2001], Wang 
et al. [2004] and DiCarlo [2004]. 

Various mathematical models have been developed to capture fingering phenomena 
[Philip., 1975; Parlange and Hill, 1976; DiCarlo et al., 2008] with many based on applying 
the Richards equation (RE) in combination with appropriate constitutive equations for 
soil properties. Techniques of linear stability analysis were applied to two-dimensional 
models by Saff'man and Taylor [1958], Chuoke et al. [1959], and Parlange and Hill [1976], 
who derived stability criteria and analytical predictions for quantities such as finger width 
and velocity. Raats [1973] postulated a criterion for stability which stated that a wetting 
front is unstable if the velocity of the front increases with depth; this is clearly satisfied for 
some layered media as well as for water-repellent soils. Many analytical results have been 
compared to experiments by authors such as Glass et al. [1989b] and Wang et al. [1998], 
who found that no single analytical formula is capable of capturing the behaviour of the 
majority of soils. Other modifications and improvements to the theory have appeared more 
recently, such as Wang et al. [1998] who modified the work of Parlange and Hill [1976] 
to include dependence on the water- and air-entry pressures. A comprehensive review of 
stability results, including comparison to experiments, can found in de Rooij [2000]. There 
has also been a great deal of recent effort on explaining gravity- driven fingering using 
models based on conservation laws [Eliassi and Glass, 2002; Nieber et al., 2003, 2005; 
Gueto-Felgueroso and Juanes, 2008, 2009a]. A recent paper by Gueto-Felgueroso and 
Juanes [2009b] presents the first exhaustive stability analysis of a conservation law that 
leads to fingering in unsaturated fiow, and a follow-up study by the same authors performs 
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an extensive comparison to experiments as well as providing an excellent review of the 
current literature [Cueto-Felgueroso and Juanes, 2009a]. 

There has been a recent surge of interest in modelling fingering instabilities using ex- 
tensions of the RE model, such as the work of Cuesta et al. [2000] and Cuesta and Hulshof 
[2003] who analyze non-monotonic travelling wave profiles that arise when dynamic cap- 
illary effects are incorporated. Both Eliassi and Glass [2001] and Nieber et al. [2005] 
identified a number of mechanisms that could give rise to gravity-driven fingering, includ- 
ing non-monotonicity in hydraulic properties, dynamic capillary effects and hysteresis. 
Egorov et al. [2003] provide an overview of the mathematical formulation showing that 
Richards equation is unconditionally stable even for heterogeneous media. Furthermore, 
Nieber et al. [2003] claim that dynamic (or non-equilibrium) effects are sufficient to cause 
formation of fingered flow, while persistence of fingers is dominated by hysteresis. In par- 
allel with these developments, several novel mathematical models have been developed 
which incorporate these and other effects. A number of authors have investigated the use 
of non-equilibrium effects [Mitkov et al., 1998; Cuesta et al, 2000; Hassanizadeh et al., 
2002; Helmzg et al, 2007; Manthey et al, 2008] while others [DzCarlo et al, 2008] have 
been inspired by non-monotonicity to introduce extra terms in the RE that capture the 
"hold-back-pile-up" effect examined by Eliassi and Glass [2003]. Non-equilibrium effects 
have also been studied in the context of two-phase flow by van Duijn et al. [2007], who 
used an extension of the Buckley-Leverett model to obtain non-monotonic profiles with 
both infiltration and drainage fronts. Sander et al. [2008] proposed a one-dimensional 
RE model including hysteresis and non-equilibrium capillary terms, which is very closely- 
related to the model studied in this paper. Adding to the controversy are experimental 
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results such as DiCarlo [2007] which failed to find significant dynamic effects in gravity- 
driven infiltration. Notwithstanding the extensive literature on this subject, many open 
questions remain about which governing equations and constitutive relations are most 
appropriate for capturing preferential fiows. 

We will focus on a specific model called the relaxation non-equilibrium Richards equa- 
tion (or RNERE) [Nieber et ai, 2003, 2005] which incorporates both dynamic and hys- 
teretic effects. These authors developed an iterative algorithm for integrating the gov- 
erning equations numerically, and showed that their method is capable of generating 
finger-like solutions. The main drawback of this work was that it contained no concrete 
comparisons to experimental results. In this paper, we perform a more extensive suite 
of numerical simulations with the RNERE model and compare the results to published 
experimental studies. We also carry out a careful numerical convergence study and in- 
vestigate the sensitivity of the model to changes in physical parameters and algorithmic 
aspects such as the choice of inter-block averaging for hydraulic conductivity. The results 
demonstrate that the RNERE model is capable of reproducing realistic fingering fiows for 
a wide range of physically relevant parameters. 

2. Mathematical Model 

We begin by presenting the governing equations for the RNERE model as presented 
by Nieber et al. [2003], while at the same time reviewing earlier work on fingered fiow in 
porous media. The RE is written in mixed form as 

^ = v*-(r(r)W)-^, (1) 
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where t* represents time [s], 6* is the volumetric water content or saturation [m^/m^], and 
ip* is the water pressure head [m]. The asterisks are used here to indicate dimensional 
quantities, and will be dropped shortly when the equations are non-dimensionalized. Hy- 
draulic conductivity [m/s] is denoted by k*{6*), which is assumed to be a given function 
of water content in the case of unsaturated flow. The spatial domain is two-dimensional 
with coordinates {x*,z*), where z* represents the vertical direction and is measured pos- 
itive downwards and x* is measured horizontally. This form of the RE is called "mixed" 
because both saturation and pressure appear as dependent variables, and it is preferred 
to both the ^-based form (which becomes singular when the flow is fully saturated) and 
the i/j-hased form (which leads to large mass conservation errors when discretized) [Celia 
et ai, 1990]. 

Many models of flow in porous media combine the RE with an equilibrium constitutive 
relation of the form ip* = p*{6*), which assumes that transport properties relax instan- 
taneously to their equilibrium values as water content varies during a wetting or drying 
process. This is a reasonable approximation under certain circumstances; however, there 
is now evidence from both laboratory experiments [DiCarlo, 2004] and stability analyses 
[Nieber et ai, 2005] that suggests the RE by itself is unable to capture the non-monotonic 
proflles observed in flngering instabilities and so it lacks some critical physical mecha- 
nism. An illustration of a typical solution proflle is shown in Fig. 1, wherein a downward- 
propagating flnger is led by a nearly-saturated "tip" region that leaves behind it a "tail" 
region having a lower water content. An earlier attempt at simulating fingered flow using 
the equilibrium RE was made by Nieber [1996] who incorporated hysteretic effects and 
found that flngers only appeared when a downwind-weighted mean was used for hydraulic 
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conductivity. Eliassi and Glass [2001] concluded that the finger-hke profiles obtained in 
these simulations did not represent solutions of the actual model equations, but rather 
were numerical artifacts arising from local truncation errors due to the particular choice 
of downwind mean. 

Based on physical arguments, Hassanizadeh and Gray [1993] advocated that a non- 
equilibrium version of the capillary pressure relationship should be employed in situations 
where the relaxation time is comparable to other time scales in the flow. This work inspired 
Nieber et al. [2003] to propose their RNERE model in which the RE was supplemented 
by a relaxation equation of the form 

* an* 

r-p*{e*) = -^, (2) 

pg ot* 
where p*{6*) represents the equilibrium water pressure head [m], p is the density of water 
[kg/m^], g is the gravitational acceleration [m/s'^], and r* = t*{iIj*,6*) > is a suitably- 
chosen capillary relaxation function [kg/ms]. They presented numerical simulations of 
finger-like instabilities and concluded that a non-equilibrium effect is sufficient to initiative 
the instabilities and that hysteresis is necessary to sustain the fingers once formed. In 
most cases, r* is assumed either to be a constant [Hassanizadeh et al., 2002; Manthey 
et al., 2008] or else a separable function of dynamic capillary pressure and water content 
[van Duijn et al., 2004; DiGarlo, 2005; Nieber et al., 2005], although the proper choice 
of functional form for the relaxation function remains an open question. We note that 
Eq. (2) should be viewed as an equation for the dynamic capillary pressure ip* rather than 
an evolution equation for 6*; indeed, when (2) is substituted into Eq. (1), the resulting 
PDE takes the form of a third-order evolution equation for 6* which is of pseudo-parabolic 
type [King and Guesta, 2006]. 
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Before proceeding further, we briefly mention several other attempts at incorporating 
non-equilibrium effects into the RE in more general contexts not directly related to fin- 
gering. Mitkov et al. [1998] took a phase-field model for solidification and adapted it to 
porous media flow; their model contains a phenomenological term in which the constants 
have no direct relationship to the physics. Barenhlatt et al. [2003] suggested an alternate 
approach in which dynamic effects are incorporated into both capillary pressure and hy- 
draulic conductivity through an "effective saturation" variable. Three further variants of 
the RE called the hypo-diffusive, hyperbolic and mixed forms were proposed with an aim 
to reproducing the "hold-back-pile-up" effect observed in experiments [Eliassi and Glass, 
2003; DiCarlo et al., 2008]. Analytical and numerical results suggest that many of these 
approaches show promise, but the proper choice of model remains an open question. 

To complete the mathematical description of the RNERE model equations (1) and (2), 
the equilibrium pressure p* and hydraulic conductivity k* must be specified as functions 
of water content. These quantities are customarily expressed in terms of the normalized 
water content 

fs — t>r 

where 6s and 6r are the saturated and residual (or irreducible) water contents, respectively; 
6 is commonly referred to as the effective saturation or simply saturation. In a partially 
saturated porous medium, the saturation variable satisfies ^ 6r ^ 6* ^ 6s ^ 4>, where 
represents the porosity, so that 6 always lies between and 1. We adopt the widely-used 
van Genuchten-Mualem relationships for p*{6) and k*{6) [van Genuchten, 1980], which 
are monotonic functions containing several empirical fitting parameters that are used to 
fit with experimental data for a variety of soil and rock types. Saturation and pressure 
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are related at equilibrium by 

d = {l + a}\pTT'"', (4) 

where n^ and m^ = 1 — 1/n^ are parameters that govern the shape of the capillary curves. 
In practice, ^ is a hysteretic function wherein the inverse capillary length a*^ [m~^] and 
shape parameter ni differ depending on whether the current state is evolving along the 
main wetting curve (£ = w) or main drying curve {i = d). The corresponding hydraulic 
conductivity is given by 

k*{9) = koVe [1 - (1 - ^1/™^)"^^] ' , (5) 

where ko represents the fully saturated value (at ^ = 1). The RE (1) becomes degenerate 
when k*{0) = 0, although we never have to deal explicitly with this issue because we always 
choose a value of initial saturation 6i greater than the residual value 6r- Furthermore, 
we have not encountered values of the reduced saturation in excess of 1 (for which the 
conductivity is undefined in Eq. (5)) and so the computed solution remains within the 
physical limits 0^9^ 1. Details of the precise form of the hysteresis model to be used 
will be provided later in Section 3. 

2.1. Non-Dimensionalization and Choice of r 

In this section, the equations are reduced to dimensionless form using the transforma- 
tions 

i) = al,ij\ p = alp*, k = k*/ko, (6) 

t = alkX/iO,-9r), 
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where (a^)^^ (the reciprocal of the van Genuchten parameter for the main wetting curve) 
has been used as the natural length scale. The governing equations (1), (2), (4) and (5) 
then reduce to 

39 

ij = p + T{^,e)—, (8) 

9 = Seip):={l + ae\pn-"'', (9) 

k{e) = Kg{e) ■.= Ve{i-{i- ^v™.)™^)^ , (lo) 

where the dimensionless capillary relaxation function is 

P9 
and S£{p) and Kf{9) represent the hysteretic constitutive relations (which depend on the 
wetting/drying state i = w,d). It will prove convenient when describing the numerical 
algorithm to recast the time-derivative in Eq. (8) in terms of the equilibrium capillary 
pressure as 

ri^,e)^ = ^-p, (11) 

where r(V^,^) = T{ip,e) (de/dp). 

Following Nieber et al. [2005] (who was motivated by the experiments of Selker et al. 
[1992b]) we assume that r is a function of ip only 

r(z^) = r„[^-z^JX, (12) 

where Tq, ipo and 7 are constants and [■]+ = max(-, 0). Nonetheless, the appropriate choice 
of functional dependence for r on the state variables remains an open question. Various 
other functional forms have been proposed by DiCarlo [2005] and Sander et al. [2008], all 
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of which we find leads to similar fingering patterns provided that r —)> as ■?/'—)■ ■i/'o and 
that the magnitude of the relaxation parameter is comparable. On the other hand, if r is 
taken to be a constant then no fingers were observed in our simulations, hence suggesting 
that it is essential to have a solution-dependent r. 

2.2. Boundary and Initial Conditions 

Inspired by the geometry most commonly employed in experimental studies, we consider 
a two-dimensional rectangular domain as shown in Fig. 2 that has width L and height 
H (both dimensions having been non-dimensionalized by scaling with a*^ like the other 
lengths in Eq. (6)). The initial saturation is assumed constant throughout the domain, 
d{x, z, 0) = di- No-fiux conditions are imposed along side boundaries and specified infiow 
(outfiow) conditions are given along the top (bottom) boundaries, both of which we express 
in terms of the dimensionless fiux variable 

q = {qi.q2) = -k{e)V{i^-z), (13) 

where gi and g2 represent the components of the fiux vector. This last equation is a 
statement of Darcy's law for unsaturated fiow and is rescaled according io q = q*/ko, 
where the original physical fiux q* has units of m/s. The fiuxes on the left, right, and 
bottom boundaries are given respectively by 

gi(0,z,t) = 0, qi{L,z,t) = 0, q2ix, H,t) = Qi, (14) 

where qi represents a constant gravity-driven background fiux that corresponds to the 
initial saturation 6i. The specification of such a background fiux is essential when the 
porous medium is not air-dried, such as is the case for most naturally occurring soils. In 
order to drive the formation of fingers, the fiux along the top boundary is specified by 
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some constant background flux Qi plus an infiltration fiux q that is imposed along a strip 
of width ^ di ^ L, 

t^ Q ^)= h^ + qix), if \2x -L\^ d^, ,^^. 

^^ ' ' '^ \ Qi, otherwise. ^ ^ 

Following Nieber et al. [2003], this infiltration flux is written as the sum of an average 
value Qs plus a small sinusoidal perturbation 

q{x) = qs + QsV cos i—{2x - L + di)\ , (16) 

where t] represents the amplitude of the perturbation and 27r f/di the frequency (for / 
a positive integer). We note that the size and number of flngers actually observed in 
simulations is relatively insensitive to the choice of perturbation parameters 77 and /. 

We close with a brief mention of a common result on stability of gravity-driven vertical 
inflltration flow, wherein viscous forces tend to stabilize the flow while gravitational forces 
are the destabilizing influence. Extensive work on stability has been reported by many 
authors, including Philip [1975], Parlange and Hill [1976], Wang et al. [1998] and de Rooij 
[2000]. It is known that unstable flow will occur if the hydraulic conductivity increases 
with depth, which translates into a requirement that the inflow at the top boundary 
satisfles < gj + g^ < 1. 

3. Solution Algorithm 

We next describe the algorithm developed by Nieber et al. [2003] for solving the RNERE 
problem in Eqs. (7) and (9)-(ll), which is an iterative strategy that employs a flnite 
volume spatial discretization in space and a semi-implicit time-stepping scheme. The 
domain is divided into an N^ x N^ rectangular grid, with cell dimensions Ax = L/N^ 
and Az = H/N^ in the x- and ^-directions respectively. The discrete saturation 9ij 
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approximates the solution at cell centers {{i — l/2)Ax, (j — l/2)Az), and similarly for 
the pressure head ipij. Employing an implicit backward Euler discretization for the time 
derivative in Eqs. (7) and (11) and centered second-order differences in space, the discrete 
equations become 



At 



A^ y"'^'/''^ A^ '''-'/''^ A^ ) + 

1 f, iJi,j+i - i^ij , iJij-iJij-i 

'^i,j+l/2 T— «^JJ-l/2 



z \ ' Az ' Az 



Kj+l/2 - Kj~i/2 /-.^N 

Az ' ^ ' 



and 



^i^i,]) ''\^ ''^ = i'^,3 - Pi,j, (18) 

for i = 1,2, . . . , Nx and j = 1,2,..., A^^. The time step is denoted by At and the "hat" 
notation 6 and p refers to a solution value at the previous time step. 

It is worth mentioning that although an upwind difference might normally be advocated 
for the convective (gravitational) term in Eq. (17), we have chosen to use a centered 
difference for consistency. In practical computations, we observe no difference between 
an upwind or centered difference treatment of the dk/dz term because capillary effects 
dominate in this problem. 

In contrast with the cell-centered saturation and pressure head unknowns, the hydraulic 
conductivity values ki±i/2j and kij±i/2 are located at cell edges. Since the conductivity 
depends on saturation which is not available at cell edges, it must be approximated using 
some weighted mean of nearby values of saturation; the specific choice of averaging method 
will be considered in detail in Section 4.1. The capillary relaxation function in Eq. (12) 
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is replaced by the regularized function Ts{ip) = To max ((■?/' — ■ipo)'^,^), where the cut-off 
parameter < 5 ^ 1 prevents r from becoming zero and hence avoids a singularity in 
Eq. (18). 

The difference stencils in Eq. (17) involve values of pressure head ipij for i = 0, N^ + 1 
and j = 0,Ny + 1, located at points lying one-half grid cell outside the physical domain. 
These "fictitious values" are eliminated using the flux boundary conditions (14)-(16) as 
follows. First, the flux is discretized along cell edges, with the gi-component (on side 
boundaries i = 0, A^^^) being approximated by 

Ql;i+l/2,j - -Ki+i/2,j T— , (19j 



Ax 
and the g2-component (along horizontal boundaries j = 0, A^^) by 

g2;i,i+l/2 = -ki,]+l/2 ( ^^ 11 . (20) 

Then the boundary conditions for gi and g2 are used to express fictitious point values 
in terms of known values of pressure head at interior points, which can then be used in 
Eq. (17). 

We now describe the iterative scheme for solving the nonlinear system (17), which can 
be written more succinctly in matrix-vector form as 

An + ^=0, (21) 

where and 11 are vectors containing the discrete approximations of 6ij and {ip — z)ij 
respectively, and A = A(0) is a symmetric pentadiagonal matrix whose entries are non- 
linear functions of saturation. Nieber et al. [2003] did not base their iterative solution 
strategy directly on Eq. (21) because the matrix A is not positive definite; instead, they 
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proposed the following modified iteration 

(A'^+i + D'^+^)n'^+i = D'^+^n'^ - ^^ ~® , (22) 

where v represents the iteration number and D is a diagonal matrix whose entries are 
given by 

^_ 1 dQ _S\P) d U^t + rmp \ 

Atd^ At d^ y r(^) + At J' ^ ^ 

This approach has the advantage that the iteration matrix (A + D) is both symmetric 

and positive definite and thus has much better convergence properties. We now outline 

the iterative procedure within each time step, assuming that each iteration begins with 

z/ = 0, e° = 0, ^° = § and P" = P: 

Step 1. Solve the relaxation equation 

pv+l _ p 

^ ^ At 
for P^+K 

Step 2. Update the saturation using 0*^+^ = S{P'^^^). 

Step 3. Evaluate the matrices A and D at = 0''+\ P = P""^^ and ^ = ^''. Then 
solve the linear system (22) for 11''+^ and let ^''+^ = 11'"+^ + z. 

Step 4- If the current solution satisfies the convergence criterion HII'^''"-'^ — n'^||2/||n'^||2 < 
10^®, then stop. Otherwise, increment z/ and return to Step 1. 

Following Eliassi and Glass [2001], we employ a variable time step which is initialized 
to At = 10^"^ and then replaced at the end of each time step by min(1.05At, Atmax), 
where the maximum allowable step is given by the CFL-like condition Atma^; = 
0.1 min(Aa;, A2;)/gs. This approach minimizes initial start-up errors by taking a rela- 
tively small time step initially, which then increases gradually to Atmax- The algorithm 
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just described is implemented in Matlab and uses the built-in preconditioned conjugate 
gradient solver peg to invert the linear system in Step 3. 

We stress that this is only one possible choice of algorithm and that many other strate- 
gies have been proposed for solving the coupled system of equations for saturation and 
capillary pressure. For example, Cuesta [2003] and Cuesta and Pop [2009] have analyzed 
a number of algorithms (including the one described above) in the context of the Burgers 
equation, supplemented by dynamic capillary effects. 

An important aspect of our RNERE model is the singularities that occur when 9 = 
(where the iteration matrix (A + D) fails to be positive definite) and 6 = 1 (where 
the derivative of the hydraulic conductivity k'{6) becomes unbounded). A number of 
methods have been proposed in the literature (e.g., Starke [2000]; Pop [2002]) to regularize 
coefficients in the governing equations in order to avoid these singularities. We have not 
made use of any such regularization here because in practice, we find that the computed 
saturation never reaches the limiting values of or 1. Nonetheless, it may be worthwhile in 
future to consider implementing such a regularization approach to improve the efficiency 
and robustness of the algorithm in cases where conditions approach the saturated and 
unsaturated limits. 

3.1. Implementation of Hysteresis 

An integral component of the RNERE algorithm is the specification of the hysteretic 
state, which is potentially different at every point in the domain and depends on the 
local saturation and wetting history. We have chosen to implement a closed-loop hystere- 
sis model described by Scott et al. [1983] and implemented by Eliassi and Glass [2003] 
wherein all curves have the same values of residual and saturated water content (0 and 
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1 respectively in our dimensionless variables). We have also taken the shape parameters 
for the wetting and drying curves to be constant and equal (n := n^ = n^), so that the 
various curves differ only in their value of a^ (although in general, the values of n^ should 
also depend on the current hysteretic state). 

Following Eliassi and Glass [2003], the main drying curve is written 6 = Sd{p) while 
the scanning drying curves are given by the scaled equation 

= ^^, (24) 

where 6 and p denote respectively the saturation and pressure reversal points along the 
previous wetting curve. Similarly, the main wetting curve is ^ = Sw{p) while the scanning 
wetting curves are written 

e = Orev + (1 - erev)SM, (25) 

where 

and 6 and p are the reversal points along the previous drying curve. 

The main drying and wetting curves are unique, while the scanning curves differ de- 
pending on the reversal points which are determined as follows. At each point in space, 
we maintain the current state (wet or dry) as well as the previous reversal point {0,p). To 
avoid problems with convergence in the iterative scheme, the hysteretic state is updated 
only at the end of each time step and not within a i/-iteration. To detect a reversal point 
along a local drying or wetting curve, we test whether the time rate of change of satura- 
tion has reversed sign between the current (k) and previous {k — 1) time steps, which is 
equivalent to checking that '^O'^j ■ ^0^^^ < where '^O'^j = 9^^ — Ofj^. To avoid spurious 
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wet/dry oscillations between successive time steps, we impose the additional constraint 
that |A^f I > £, where e is a reversal threshold. If both of these criteria are met, then a 
flow reversal has occurred and the current state is switched to either wetting (if A6^j ^ 0) 
or drying (if A^f • < 0), and the current values for the reversal point {0,p) are updated. 
The appropriate scanning curve - either Eq. (24) or (25)-(26) - is then used to determine 
the capillary pressure as a function of saturation. 

4. Numerical Simulations 

To investigate the relevance of the proposed model and the accuracy and efficiency of 
the numerical algorithm, we consider a "base case" corresponding to a 14/20 grade sand 
studied experimentally by Glass et al. [1989b]. The parameters listed in Table 1 are taken 
directly from their paper, with the exception of n, 6i and a^, whose values are justified 
in Sections 4.4-4.6. 

The computational domain is taken to be a rectangle of width L = 14 and height H = 35 
in dimensionless units, where L is chosen slightly larger than the actual infiltration width 
di = 10.5 used in experiments in order to minimize boundary effects. Unless otherwise 
noted, the domain is discretized using a uniform grid having A^^: = 201 and Nz = 401 
points in the x and z directions respectively. All simulations were performed on a Mac 
Pro with 2x3 GHz processor and 8GB RAM, with a typical run requiring approximately 
2 hours of computation time. 

A sample computation with the base case parameters is shown in Fig. 3, which depicts 
the progression of the wetting front at a sequence of equally-spaced times. In these 
plots (as well as the other plots and tables that follow) all quantities are expressed in 
dimensionless form. The plotted contours of saturation correspond to a value of 6 equal 
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to 25% of the finger tail saturation, whicli we liave found gives a good representation of 
the finger size and shape. The structure of the individual fingers is seen more clearly 
in the saturation map given in Fig. 4, where each concentrated finger tip is followed 
by a tail region of roughly constant saturation. The "capillary fringe" region, depicted 
schematically in Fig. 1, is evident as a narrow zone of rapid saturation change surrounding 
each finger. The finger tip and tail are evident in this plot and the shape of each finger 
is in qualitative agreement with the generic profile sketched in Fig. 1. 

To illustrate the importance of dynamic and hysteretic effects in the RNERE model, 
we show in Fig. 5 how the wetting front differs when either of these effects is left out. 
When the dynamic term is omitted from the saturation equation (see Fig. 5(a)) fingering 
instabilities clearly fail to be initiated. Conversely, when hysteretic effects are left out 
(see Fig. 5(c)) protrusions begin to form at the wetting front but they never actually 
develop into full-blown fingers. These observations are consistent with the claim of Nieber 
et al. [2003] that dynamic capillary effects are responsible for the initiation of fingering 
instabilities, while hysteresis is required to sustain the fingers in time. 

In the following sections, we present an extensive suite of numerical simulations that 
address the following four issues: 

• choosing an appropriate mean for the estimation of inter-block hydraulic conductivity 
values; 

• determining the dependence of the numerical solution on grid resolution, and com- 
paring to previously published simulations; 

• measuring the sensitivity of the solution to certain key parameters: shape parameter 
(n^), dynamic relaxation coefficient (r) and initial saturation (^j); and 
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• comparing simulated results to previously published experimental data. 
In addition to providing plots of saturation plots, we will also report quantities such as 
finger width (df), finger velocity (vf), number of fingers (Nj), and average volume fiow 
rate through each finger (Qf), all of which vary depending on the value of the infiltration 
fiux Qg. When multiple fingers are present and a specific quantity varies from finger to 
finger, we report the average value over all fully-developed fingers. 

4.1. Choice of Mean for Inter-Block Conductivity 

As mentioned in Section 3, the discrete equations require values of hydraulic conductiv- 
ity at cell edges (/ci±i/2,j and kij±i/2) whereas the values of saturation on which k depends 
are defined at cell centers; therefore, some form of averaging is usually necessary. It is 
well known that discrete approximations of the RE can be very sensitive to the choice 
of inter-block averaging used for hydraulic conductivity [Belfort and Lehmann, 2005]. A 
number of different approaches have been advocated in the literature, for instance us- 
ing arithmetic [van Dam and Feddes, 2000], geometric [Haverkamp and Vauclin, 1979], 
harmonic [Das et al, 1994], and Darcian- weighted means [Warrick, 1991]. Cardwell and 
Parsons [1945] showed that the effective permeability for a heterogeneous porous medium 
must lie somewhere between the harmonic and arithmetic mean values. Furthermore, 
Warren and Price [1961] used Monte Carlo simulations of random media to show that 
the expected value of conductivity for a heterogeneous system is given by the geometric 
mean. 

Of particular interest in this paper is the case where the conductivity undergoes large 
variation between grid cells owing to the presence of sharp wetting fronts at finger bound- 
aries. Although a straightforward analytical argument indicates that the harmonic mean 
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is the appropriate mean to use in such situations in ID [Gutjahr et al, 1978], this is 
not the case in higher dimensions where many computational studies indicate that the 
harmonic average yields results that are inferior to those using other averaging methods 
[Haverkamp and Vauclin, 1979; Belfort and Lehmann, 2005; Finales et al., 2005]. Ex- 
tensive comparisons have been drawn using measures such as resolution and stability of 
the wetting fronts, sensitivity to grid refinement, and robustness over a wide range of soil 
types and physical parameters. There remains a significant degree of controversy over 
which averaging procedure is best in practice, and to date no single mean has been found 
to be superior in all circumstances. 

In this section, we compare results using the arithmetic and geometric means for hy- 
draulic conductivity, which we have found are the most common means utilized in com- 
putations. Values of conductivity along vertical cell edges are determined as follows 

Arithmetic mean: ki±i/2j = -{ki±ij + hj), 
Geometric mean: ki±i/2j = \/ki±ijkij , 

with similar formulas for kij±i/2 along horizontal cell edges. Problem parameters are 
taken from simulations presented by Nieber et al. [2003], which are identical to those for 
our base case described earlier, except for a few differences indicated in Table 2(b). The 
same set of parameters will also be considered in the following two sections. Our model is 
identical to Nieber et a/.'s, except for a slight difference in the implementation of capillary 
hysteresis. 

The results for the geometric mean are shown in Fig. 6 for five different choices of grid 
resolution (51 x 61, 101 x 121, 201 x 241, 401 x 481 and 801 x 961), and the solution has 
clearly converged to a wetting profile with 4 well-defined fingers even on a 201 x 241 grid. 
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In contrast, Fig. 7 demonstrates that simulations with the arithmetic mean converge more 
slowly, and the results on the finest grid are still not fully converged. These simulations 
are consistent with those of Zaidel and Russo [1992], who found that use of the arithmetic 
mean introduces excessive smearing in wetting fronts and underestimates saturation values 
relative to the geometric mean. Based on these results, we conclude that the geometric 
mean is superior for the problem under consideration, which is also consistent with the a 
number of previous studies [Haverkamp and Vauclin, 1979; Hornung and Messing^ 1983; 
Belfort and Lehmann, 2005]. Consequently, we have chosen to apply the geometric mean 
in all remaining computations in this paper. 

4.2. Grid Refinement Study 

To ensure that the numerical solution does converge with the expected second order 
accuracy, simulations were performed on a sequence of successively refined grids of size 
101 X 121 to 801 X 961 cells. The latter represents the finest resolution possible owing 
to memory restrictions on the computing equipment readily available to us. The grid 
resolution and physical parameters in this case were chosen to correspond to the numerical 
simulations of Nieber et al. [2003]. 

The solution on the finest grid is treated as the "exact solution" and the error is esti- 
mated using the £2 norm of the difference between exact and computed values of satura- 
tion. The resulting absolute errors are summarized in Table 3 from which it is clear that 
the solution converges as the grid is refined; furthermore, the order of convergence is close 
to the expected value of 2. Fig. 8 depicts saturation contours corresponding to various 
grid refinement levels and clearly demonstrates the convergence of the numerical solution. 
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4.3. Comparison with Nieber et al.'s Computations 

We now investigate the effects of grid resolution in more detail by way of a direct 
comparison with simulation results reported by Nieber et al. [2003]. We focus on their 
second set of simulations (c.f., their Fig. 8), using different values of the infiltration width 
{di = 1, 5, 10.5, 15, 20, 25 and 30) and grid resolutions of 101 x 121 and 201 x 241. Our 
results are summarized in Fig. 9 from which we observe a close match with Nieber et al. 
[2003, Fig. 8] in terms of both number of fingers and saturation levels. There are slight 
differences between finger widths and velocities, but we attribute these to our alternative 
implementation of capillary hysteresis. We emphasize the discrepancies between our high 
and low resolution results for values of di = 10.5, 15 and 30, which stresses the importance 
of using a sufficiently resolved grid in these fingering computations. It is particularly 
important to use a higher grid resolution in situations such as di = 30 that are close to a 
"transitional phase" where in this case the solution exhibits somewhere between three and 
four fingers. We also stress the importance of performing a careful convergence analysis 
as part of any numerical study of fingering to ensure that the features being simulated 
are a true representation of actual fingering instabilities of the governing equations. 

4.4. Sensitivity to Capillary Shape Parameter, n 

In the next three sections, we switch to the base case and investigate the sensitivity of the 
solution to changes in a number of important parameters. No value is provided by Glass 
et al. [1989b] for the capillary shape parameter n appearing in the van Genuchten-Mualem 
relationships (9) and (10), and so we look for guidance in related experimental studies 
on sandy soils. The values for n reported in the literature exhibit significant variability, 
lying anywhere between 3 and 20 even for porous media having similar coarseness and 
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wettability (e.g., Schroth et al. [1996]; Nieber et al. [2000]). The larger values of n typically 
correspond to water-repellent soils in which fingers are more likely to form, while smaller 
values indicate a reduced tendency to generate fingering instabilities. Consequently, it is 
important in any modelling study of fingering to understand the effect of changes in n on 
the character of the solution. 

In Fig. 10, we present simulations for several values of n lying between 4 and 15, while 
all other parameters are set to the base case values. In all simulations except for n = 4, the 
solution exhibits well-defined fingers having the characteristic non-monotonic saturation 
profile down the central axis of each finger. Furthermore, increasing n leads to an increase 
in the tip/tail saturation ratio and a decrease in finger width and finger fiux, as expected; 
in other words, fingers are more diffuse for smaller n while finger boundaries become 
sharper when n is increased. The value of n = 12 used in the base case is characteristic of 
highly water-repellent soils, which is corresponds to the other experimental and numerical 
studies we are most interested in. 

4.5. Sensitivity to Initial Saturation, 6i 

A study of the effect of the initial saturation 9i on our numerical solution is warranted 
for two reasons. First of all, the nature of fingering instabilities and the properties of 
individual fingers (such as finger width and velocity) can be very sensitive to the choice of 
initial water content, as evidenced by several experimental studies [Diment and Watson, 
1985; Banters et al., 2000; Wang et al., 2003]. Finger shape and size depends strongly 
on the initial wetting state; in particular, vertical infiltration into a soil with larger ini- 
tial saturation tends to generate fingers that are more diffuse than when the soil is dry. 
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Secondly, as with the capillary shape parameter n, the value of 9i is frequently omitted 
in the list of parameters reported in experimental studies (e.g., Glass et al. [1989b]). 

We therefore perform a series of simulations using various choices of initial saturation 
between 0.001 and 0.075, holding the infiltration flux g^ constant. The results are sum- 
marized in Fig. 11 from which we observe that the number of fingers increases as 9i is 
increased. In fact, the spacing between fingers also decreases to the extent that when 
di ^ 0.05, the individual fingers merge together to form a single finger. As 9i and finger 
width increase, we notice from the saturation maps in Fig. 12 that the maximum fin- 
ger tip saturation decreases while the finger velocity remains relatively unchanged; this 
behaviour can be justified using a simple mass conservation argument. Most of these 
computed trends are consistent with experiments, the exception being the finger velocity 
for which some experimental studies exhibit a stronger dependence on 9i (e.g.. Banters 
et al. [2000]). 

We mention in closing this section that in the absence of a given value of initial satu- 
ration in Glass et al. [1989b], we have chosen 9i = 0.01 for the base case. This value lies 
within with the typical range of residual saturations seen in experiments for similar soils, 
and also generates fingers with a tip saturation that is consistent with values reported by 
DiGarlo [2004]. 

4.6. Sensitivity to Capillary Relaxation Coefficient, To 

Although several recent models for gravity-driven fingering have used a capillary relax- 
ation term to incorporate dynamic effects, there remains a great deal of uncertainty in 
both the functional form and overall magnitude of the relaxation coefficient r [Juanes, 
2008; Manthey et al., 2008]. Stauffer [1978] derived an empirical estimate based on the 
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Brooks-Corey model for conductivity and capillary pressure (in lieu of the van Genucliten- 
Mualem relationships used here) that takes the following form 

" = ^^' ^''^ 

where 7^ = 0.1 is a fitting parameter, /i is the fluid viscosity, and A and hb are the Brooks- 
Corey parameters. Recent experimental results suggest that the values of r for sandy 
media can range between 0.006 and 20, so that r* lies between 2 x 10^ and 6 x 10^ kg/m s 
[Manthey et al, 2008]. 

We have run a number of simulations using the functional form for r given in (11) and 
with the scaling constant To varying between 0.01 and 1.0. The resulting saturation is 
depicted in Fig. 13, from which it is evident that To has a strong influence not only on 
the number of fingers but also on finger width and velocity. If r^ is taken very small (less 
than 0.001) then dynamic effects become negligible and finger formation is suppressed. 
If, on the other hand. To is taken larger then the finger tip saturation tends to increase 
which in turn reduces the number of fingers. 

For the purposes of the base case, we have chosen an intermediate value of Tq = 0.1 
which gives a range of r that is centered on the empirical estimate in Eq. (27), and which 
also corresponds well to the range of experimental values reported in the literature. 

4.7. Comparison with DiCarlo's Experiments 

In this section we consider the experimental results reported by DiCarlo [2004], who 
studied the importance of non-equilibrium effects on finger formation in sandy porous 
media. These experiments investigated the effect of changes in infiltration fiux, initial 
saturation, and porous media properties on the resulting saturation profiles. DiCarlo 
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also proposed an RE-based model which neglected hysteretic effects, but did include 
a dynamic capillary term (as in our Eqs. (7)-(8)) with a number of different forms 
for the dynamic relaxation coefficient t{6), including a constant and various power-law 
forms similar to Eq. (12). The correspondence between his numerical simulations and 
experiments (in terms of finger tip saturation) was less than satisfactory; in particular, 
although a reasonable fit was obtained for the tip saturations when r was a power law 
function, the wetting front ahead of the finger tip was much too diffuse. 

We focus primarily on DiCarlo^s experimental and numerical results for 20/30 sands 
which are reproduced in Fig. 14(a). These results were essentially one-dimensional be- 
cause the diameter of the soil columns being studied was less than the characteristic finger 
width and hence was too small for fingers to form; we have therefore performed a "quasi- 
ID" simulation in which the horizontal extent of the domain is only two grid points wide. 
Here, we choose parameters the same as in the base case except that the capillary shape 
parameter, initial saturation, and domain size are modified according to Table 2(c). Sim- 
ulations were performed for a range of values of infiltration fiux g^, and the resulting tip 
and tail saturations are plotted in Fig. 14(a). 

We extracted values of all parameters from DiCarlo [2004] except for the initial satura- 
tion which he did not provide. The computed tip saturation is quite sensitive to 6i, while 
the impact on tail saturation is much less. As 6i is increased, the tip profile in Fig. 14(a) 
shifts toward the right until it eventually overlaps with the tail profile, which corresponds 
to a stable flow. On the other hand, as 6i is decreased the profile shifts to the left and 
steepens becoming similar in shape to the "DiCarlo tip" curve. Consequently, we have 
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used initial saturation as a fitting parameter and chose a value of di = 0.001 that yields 
the best match with experimental tip data. 

The results in Fig. 14(a) demonstrate a significant improvement over DiCarlo^s model, 
especially in terms of the tip saturation. There is also excellent agreement with the tail 
data, although unfortunately there is no corresponding tail simulation from DiCarlo for 
us to compare to. A second comparison is made for a 30/40 sand from DiCarlo [2004] in 
Fig. 14(b). Except for some small deviations at the lowest infiltration rates, these results 
also show a good fit between our model and DiCarlo^s experiments. Similar comparisons 
are obtained for other soil types. 

Finally, it is worth emphasizing that our computations exhibit fingers that sustain a 
sharp front ahead of the finger tip, exhibiting none of the non-physical diffusive smoothing 
observed in the model results of DiCarlo [2005] . This discrepancy can be justified following 
Nieber et al. [2003] who attributed the initial formation of fingers to dynamic capillary 
effects that are present in both DiCarlo's and our RNERE model; however, fingers persist 
in time only when hysteretic effects are also incorporated, which is the case for our RNERE 
model but not DiCarlo^ s. 

4.8. Comparison with Glass et al.'s Experiments 

We next make use of the dimensional analysis and experiments of Class et al. [1989a, b] 
to assess the response of finger width and tip velocity to changes in infiltration fiux, 
Qs. Their two-dimensional experiments involved two layers of fine-over-coarse sand, with 
water fed in from the top and air allowed to escape freely. Fingering was observed in the 
lower, coarse sand layer - a result that can be predicted using the stability analysis of 
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Raats [1973] . We therefore restrict our attention to the lower layer only which contains a 
coarse 14/20 silica sand having grain diameter in the range 0.00070-0.0012 m. 

The parameter values used in this section are the same as those for the base case listed 
in Table 2(a). The residual saturation 6r = 0.078 is consistent with the tail water content 
reported in Glass et al. [1989b]; however, they did not provide values of the remaining 
porous medium parameters and so we chose n = 12, 6i = 0.01 and a^ = 35 m~^, which 
are consistent with other 14/20 sands in the literature. 

In Fig. 16, we present saturation contours from a series of simulations in which the 
infiltration flux g^ is varied between 0.038 and 0.32 cm/min. Decreasing the g^ causes 
an increase in the number of fingers, in addition to decreasing both finger velocity and 
tip saturation (as indicated in Fig. 15). The corresponding numerical values for various 
quantities are summarized in Table 4. 

Following Glass et al. [1989a], we relate the average finger velocity (f/) and width (df) 
to the average volume flow rate in a finger using Qf = dfVf. Since infiltration flux is 
related to finger velocity via g^ = NfdfVf/di, the finger volume flow rate can be written 
as 

Qf = ^. (28) 

We then take our simulations for different values of Qs and redisplay the results in Fig. 17 
as a plot of finger velocity versus finger flow rate, including data from Glass et a/.'s 
experiments. There is very close agreement between the simulated and experimental 
results. In particular, as the infiltration flux Qs increases (or equivalently, Vf increases) 
the slope of the velocity-fiux curve decreases. 
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In Fig. 18, we present a plot of finger width versus flow rate in which the dependence 
is approximately linear. This behavior is consistent with Glass et a/.'s experiments where 
they used a linear least squares fit to predict the finger width. However, the correlation 
here is not as strong and our computations significantly over-predict the finger width at 
higher values of finger flux. In an effort to explain this discrepancy, we plot finger velocity 
against finger width in Fig. 19, which includes the experimental data of Glass et al. 
[1989b]. The experimental points are classified as corresponding to "side" and "inner" 
fingers (where side fingers lie immediately adjacent to the side boundaries) and significant 
differences are apparent between the two sets of fingers which Glass et al. attribute to 
boundary effects. If we focus only on the interior fingers, then our model does a very good 
job of capturing the observed behaviour. Indeed, it is the contribution of the side fingers 
to the average finger width that leads to the deviations in slope at higher flux in Fig. 18. 

5. Discussion 

In this study we investigated the ability of the RNERE model to capture gravity-driven 
fingering in unsaturated soils. Our results expand on previous studies of the RNERE 
model in two ways: first, by performing an extensive sensitivity analysis for various im- 
portant physical parameters; and second, by drawing a systematic comparison between 
simulations and previously published experimental results. We showed that with a care- 
ful choice of initial saturation and capillary relaxation coefficient, the model is capable 
of accurately reproducing the fingering behaviour observed in experiments. Comparisons 
with several independent experimental studies attest to the accuracy and robustness of 
the RNERE approach. 
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In contrast to the work of DiCarlo [2005], who concluded that the RNERE does not 
contain all the required physics to describe gravity- driven fingering instabilities, we have 
shown that by coupling both non-equilibrium and hysteretic effects it is possible to cap- 
ture fingering phenomena with the RNERE. Our numerical simulations demonstrate the 
importance of performing a detailed numerical convergence study in order to ensure that 
fingers have been sufficiently well resolved. The model sensitivity analysis showed that 
dynamic capillary terms must be properly handled if the fingered fiow is to be captured 
accurately, and in particular that an accurate estimate of the Tq parameter is essential. As 
more research is undertaken in the study of non-equilibrium capillary effects, we expect 
that more accurate and reliable experimentally-validated correlations for the capillary 
relaxation parameter will become available. 

There are a number of possible avenues for future work that will be explored: 
• We will investigate the use of alternate iterative strategies that improve on the ro- 
bustness and efficiency of the RNERE algorithm. We hope to draw inspiration in this 
respect from other well-known work on RE-based methods such as Celia et al. [1990] and 
Miller et al. [1998]. Current advances in ODE solvers for dealing with event detection and 
non-smooth or discontinuous coefficients may also yield improvements in the treatment 
of hysteretic switching criteria, which has a big impact on convergence of the iterative 
scheme. Furthermore, there has been an explosion of recent work on alternate models for 
handling dynamic capillary effects which could be applied here [Beliaev and Schotting, 
2001; Sander et al., 2008; Peszynska and Yi, 2008; Helmig et a/., 2007; Manthey et al, 
2008]. 
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• Analytical results for fingered flow, derived using asymptotic or other approximate 
methods, will be studied to gain a better understanding of the impact of hysteresis and 
dynamic effects on the mechanics of finger formation. We will initially be guided by other 
previous work on traveling wave approximations for wetting fronts in the RNERE model 
with dynamic capillary effects [DiCarlo et al, 2008; Nieber et al, 2005] and hysteresis 
[Dautov et al, 2002; Sander et al, 2008]. 

• Two alternate mathematical models have recently been proposed that are relevant 
to capturing gravity-driven fingering phenomena. The model of Cueto-Felgueroso and 
Juanes [2008] accounts for effective surface tension phenomena due to saturation gradi- 
ents through the addition of a new fourth order derivative term in the equations. Their 
numerical results capture the main qualitative features of fingered flow without the need 
for hysteresis, although they state that hysteretic effects are still important and that hys- 
teresis can be easily incorporated into their model [Cueto-Felgueroso and Juanes, 2009a]. 
Another approach proposed by Pop et al. [2009] introduces an additional PDE for the 
interfacial area that obviates the need for an explicit treatment of hysteresis. We intend 
to perform an extensive computational comparison of these two models using a general- 
ization of our RNERE approach, initially in ID, that should help to elucidate the relative 
advantages and disadvantages of the various approaches. 
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tail 





Figure 1. Left: A typical finger propagating into a porous medium having an initially 
uniform saturation 9,.. At the leading edge of the finger is a well-defined "core" or "tip" 
region with water content close to the saturated value Og- Behind the tip lies a "tail" region 
having a nearly constant intermediate value of water content. Right: The corresponding 
saturation profile along the central axis of the finger. 
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Figure 2. The computational domain witli widtli L and lieiglit H. Zero flux conditions 
are imposed on side boundaries and a background flux of qi on the top and bottom 
boundaries. Finger formation is driven by an additional inflltration flux q applied along 
a portion of the top boundary having width di. 
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Table 1. "Base case" parameters given in SI units. With the exception of n, di and 
a^, these values are taken from Glass et al. [1989b]. 



Symbol Description Value Units 

n Capillary shape parameter 

a^ Inverse capillary length 

9i Initial water content 

9s Saturated water content 

Or Residual water content 

ko Saturated conductivity 

K Permeability 

ipwe Water entry pressure 



12 


— 


35 


m-i 


0.01 


m? /m? 


0.42 


m? /m? 


0.075 


m^ jvp? 


0.063 


m/ s 


6.5 X 10-1° 


9 

m 


-0.023 


m 
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Table 2. Dimensionless parameter values: (a) "base case" and Glass et al. comparisons; 
(b)-(c) modifications to base values for simulations indicated. 



Symbol Description 



Value 



(a) Parameters from Glass et al. [1989b] - "base case' 



n 


Capillary shape parameter 


12.0 


ayj 


Inverse capillary length (wetting) 


1.0 


Old 


Inverse capillary length (drying) 


0.5 


Oi 


Initial water content 


0.01 


To 


Relaxation coefficient 


0.1 


7 


Relaxation exponent 


1 


^o 


Relaxation parameter 





e 


Hysteretic reversal criterion 


10-10 


5 


Relaxation cut-off 


0.04 


Qi 


Background flux 


3.3 X 10 


qs 


Infiltration flux 


0.14 


V 


Perturbation amplitude 


0.01 


f 


Perturbation frequency 


5 


di 


Infiltration source width 


10.5 


^end 


End time 


77 


H 


Domain height 


35 


L 


Domain width 


14 


(b) Modifications for Nieber et al. [2003] (Figs. 8 & 9): 




n 


Capillary shape parameter 


7.0 


0i 


Initial water content 


0.1 


To 


Relaxation coefficient 


5.0 


(Is 


Infiltration flux 


0.2 


'(■end 


End time 


96 


H 


Domain height 


60 


L 


Domain width 


30 



(c) Modifications for DiCarlo [2004, Tab. 1] for 20/30 and 30/40 sands (Fig. U): 
n Capillary shape parameter 6.23 / 10.0 

Oi Initial water content 0.001 

H Domain height 7.08 / 6.92 

L Domain width 0^ 

^ DiCarlo^?, experimental soil columns measure less than one finger width in diameter, 
DRAFT May 19, 2010, 2:48am D ^ A F T 

so that the flow is essentially one-dimensional; we perform "quasi-lD" simulations by 
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Figure 3. Saturation contours for tlie base case, plotted at four equally-spaced time 



intervals between t = and tend = 77. 




Figure 4. Saturation map for the base case at time tend = 77. 
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(b) 



Figure 5. The effect of hysteresis and dynamic effects on the wetting front, for the base 
case: (a) with hysteresis only; (b) base case, with both hysteresis and dynamic effects; (c) 
with dynamic effects only. 



51x61 cells 101xl21cells 201x241 cells 401x481cells 801x961 cells 
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Figure 6. Saturation maps corresponding to the geometric mean for five different grid 
resolutions, using parameters from Nieber et al. [2003]. 
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Figure 7. Saturation maps corresponding to the arithmetic mean for five different grid 
resolutions, using parameters from Nieber et al. [2003]. 



DRAFT 



May 19, 2010, 2:48aiii 



DRAFT 



X - 48 CHAPWANYA AND STOCKIE: SIMULATIONS OF GRAVITY-DRIVEN FINGERING 

Table 3. Grid refinement study, where the order of accuracy is estimated as the 
base-2 logarithm of the ratio of successive errors. The "exact" solution corresponds to an 
801 X 961 grid computation. 

No. of cells {Nx X A^^) £2— error Ratio Order 



60 X 51 



121 X 101 



241 X 201 



481 X 401 



10.59 3.52 1.82 



3.01 



3.82 1.93 



0.79 4.31 2.11 



0.18 






121 by 101 cells 

- - 241 by 201 cells 

481 by 401 cells 

-- 961 by 801 cells 
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Figure 8. Saturation contours corresponding to the different grid resolutions used in 
the convergence study. 
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(a)d. =1 (b)d. =5 (c)d. =10 (d)d. =15 (e) d. =20 (f) d. =25 (g) d. =30 
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Figure 9. Comparison of fingers for various values of the infiltration width di, based on 
parameters from Nieber et al. [2003]. The solid (red) saturation contours correspond to 
solutions on a 101 x 121 grid, while the broken (black) contours are for a 201 x 241 grid. 
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Figure 10. Contour plots of saturation for different values of the capillary shape 
parameter n, with contours shown at four equally-spaced times. The base case is indicated 
by a dagger (f). 



DRAFT 



May 19, 2010, 2:48aiii 



DRAFT 



X-50 



CHAPWANYA AND STOCKIE: SIMULATIONS OF GRAVITY-DRIVEN FINGERING 



(a) (b) (c) (d) 

di = 0.001 di = O.Olt e, = 0.03 6», = 0.075 
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Figure 11. Contour plots of saturation for different values of the initial saturation 6i 
shown at four equally-spaced times. The base case is indicated by a dagger (f). 
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Figure 12. Saturation maps at the final time corresponding to the same values of 6i 
depicted in Fig. 11. 
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(a) To = 0.01 (b) To = 0.05 (c)To = 0.lt (d) t^ = 0.3 
T =0.01 X =0.05 X =0.1 X =0.-1 
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Figure 13. Contours of saturation for different values of the capillary relaxation 
coefficient Tq, at four equally-spaced times. The base case is indicated by a dagger (f). 
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Figure 14. Comparison of tip and tail saturations from the RNERE model and 
experiments. Results are for two different soil types - (a) 20/30 sand (top), (b) 30/40 
sand (bottom) - and the experimental data points are extracted from DiCarlo [2004, 
Figs. 6 and 8]. The "DiCarlo tip" curve in (a) corresponds to a numerical simulation 
using a non-equilibrium model without hysteresis from DiCarlo [2005, Fig. 2]. 
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Figure 15. Computed tip and tail saturations for the Glass et al. [1989b] comparisons. 
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Figure 16. Saturation contours for various values of infiltration flux fluxs, correspond- 
ing to parameters listed in Table 2(a) for Glass et al. [1989b]. A further comparison of 
specific quantities is provided in Table 4. The base case is indicated by a dagger (f) and 
the end time for each simulation is indicated on the plot. 
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Table 4. Comparison of the finger number, width and velocity corresponding to the 
simulations in Fig. 16. The experimental data are taken from Glass et al. [1989b] (no 
data were available for the highest flux value, g^ = 0.52). 

Experimental data Numerical simulations 



Nf Qf df Vf Nf Qf df Vf 



0.012 


4 


0.03 


0.45 


0.12 


7 


0.02 


0.28 


0.08 


0.038 


4 


0.10 


0.52 


0.21 


6 


0.07 


0.34 


0.18 


0.088 


5 


0.18 


0.60 


0.29 


6 


0.15 


0.49 


0.29 


0.11 


6 


0.19 


0.61 


0.30 


6 


0.19 


0.54 


0.32 


0.14t 


4 


0.38 


0.79 


0.41 


5 


0.29 


0.70 


0.38 


0.28 


6 


0.50 


0.91 


0.46 


4 


0.74 


1.34 


0.53 


0.32 


5 


0.66 


1.08 


0.52 


3 


1.12 


1.79 


0.59 


0.52 


- 


- 


- 


- 


2 


2.73 


3.88 


0.69 
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Figure 17. Plot of finger flow rate versus velocity with experimental data (square 
points) taken from Glass et al. [1989b]. 
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Figure 18. Plot of flow rate versus finger width with experimental data (square points) 
taken from Glass et al. [1989b]. 
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Figure 19. Plot of finger width versus finger velocity, with experimental data taken 
from Glass et al. [1989b] . Square data points denote fully developed interior fingers while 
crosses denote fingers adjacent to side boundaries. 
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